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VALIDITY OF HEAVY TRAFFIC STEADY-STATE 
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We consider a single class open queueing network, also known 
as a generalized Jackson network (GJN). A classical result in heavy- 
traffic theory asserts that the sequence of normalized queue length 
processes of the GJN converge weakly to a reflected Brownian mo- 
tion (RBM) in the orthant, as the traffic intensity approaches unity. 
However, barring simple instances, it is still not known whether the 
stationary distribution of RBM provides a valid approximation for the 
steady-state of the original network. In this paper we resolve this open 
problem by proving that the re-scaled stationary distribution of the 
GJN converges to the stationary distribution of the RBM, thus val- 
idating a so-called "interchange-of-limits" for this class of networks. 
Our method of proof involves a combination of Lyapunov function 
techniques, strong approximations and tail probability bounds that 
yield tightness of the sequence of stationary distributions of the GJN. 

1. Introduction. Jackson networks are one of the most fundamental ob- 
jects in the theory of stochastic processing networks, owing to the remark- 
able simplicity of the product form stationary distribution of the queue 
lengths. Product form results have since been established for wider classes 
of open and closed queueing networks, however, these results typically hinge 
on exponential distributions being imposed on the interarrival and service 
times. Relaxing these assumptions to i.i.d. interarrival and service times 
with general distributions significantly complicates the analysis and, barring 
special cases, does not lead to closed form expressions for the steady-state 
queue lengths distribution. From an application standpoint, however, such 
assumptions typically constitute a more adequate model of "real world" sys- 
tems. Consequently, the study of open single class queueing networks with 
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Bernoulli type routing and i.i.d. interarrival and service times, hereafter 
referred to as generalized Jackson networks (GJN), has been the focus of 
vigorous academic research. 

In the absence of closed form solutions for the steady-state queue length 
distribution in GJN, efforts have mostly centered on the development of 
various bounds and approximations. Most notably, much of the recent work 
on single and multiclass queueing networks has focused on fluid models for 
stability analysis, diffusion approximations for performance analysis and the 
relationship between the two methods and objectives. Below we provide a 
brief overview of these developments in the context of GJN. 

In terms of diffusion approximations, Reiman's seminal paper [33] proved 
that, as the traffic intensity approaches one, the normalized vector of queue 
length processes converges in distribution to a Brownian motion which is 
reflected at the boundaries of the positive orthant. The normalization in 
Reiman's heavy-traffic limit theorem corresponds to the functional central 
limit theorem, that is, time is accelerated linearly and space is compressed by 
a square-root factor. The limit process, known as reflected Brownian motion 
(RBM), was first constructed in earlier work of Harrison and Reiman [18]. 
The main appeal of the RBM process is that it provides a parsimonious (two 
moment) mesoscopic approximation of the underlying discrete dynamics in 
the GJN. Reiman's theorem provides a rigorous justification for the use of 
RBM as an approximation to the original (appropriately normalized) queue 
length process in the GJN. 

To approximate the steady-state behavior of the network, it is first nec- 
essary to derive conditions under which both the GJN and RBM admit a 
stationary distribution. To this end, Sigman [35] and Down and Meyn [15] 
were the first to prove that the GJN possesses a unique stationary distri- 
bution if and only if the traffic intensity at every station is less than unity. 
A similar result was later proved by Dai [12] using a remarkable connection 
between the stochastic and fluid models of queueing networks; see also [37]. 
(Fluid limits of GJN were first characterized by Johnson [22] and Chen and 
Mandelbaum [9, 10].) Despite the fact that the stability conditions in GJN 
have such a simple form, exact computation of the stationary distribution 
is not possible except in the case of product form networks. 

It turns out that the stability condition for RBM is quite similar to 
that of the underlying GJN. In particular, it was shown by Harrison and 
Williams [19] that a (unique) stationary distribution exists for the RBM, 
roughly speaking, if and only if the "traffic intensity" at each station is less 
than unity (a precise definition will be given in the following section). For 
further details and references, the reader is referred to the survey paper by 
Williams [38]. Although the stationary distribution of an RBM cannot be 
typically computed in closed form, with the exception of special cases in- 
volving a skew-symmetry condition (see, e.g., [19] and further discussion in 
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Section 4.4), it can be solved for numerically. This computation typically 
exploits the partial differential equation, also known as the basic adjoint 
relation, that characterizes the stationary distribution. The algorithms of 
Dai and Harrison [13] and Shen et al. [36] build on this idea and describe 
different approaches to compute the stationary distribution by numerically 
solving the basic adjoint partial differential equation. 

Given the apparent solidarity between GJN and its associated RBM ap- 
proximation, as well as the fact that steady-state computations can be more 
easily performed for the latter, it has become commonplace to approximate 
the steady-state of GJN with that of the associated RBM. It is important 
to note that this approximation is not rigorously justified and, in particu- 
lar, does not follow from Reiman's [33] heavy-traffic process level approx- 
imation. To date, this correspondence has only been established in some 
specific problem instances, starting with Kingman's seminal work [27, 28] 
on the GI/GI/1 queue. Specifically, Kingman proved that the distribution 
of the diffusion-scaled stationary waiting time converges to an exponential 
distribution, which is identical to the steady-state of the RBM obtained as 
a process limit from the normalized queue length. (For extensions of this 
result that cover dependent input processes, see, e.g., [34].) In the context 
of a single station multiclass queueing network with feedback, the workload 
process, after diffusion scaling, also converges to a one-dimensional RBM. 
The steady-state mean value of this RBM is known to be equal to the ex- 
pected workload in the original model under diffusion scaling. Thus, in this 
example the validity of the RBM steady-state approximation is established 
thanks to the fact that both of the above expectations can be computed 
in closed form. Harrison [17] proved that, in a series of two single-server 
queues, the scaled steady-state vector of queue lengths converges weakly to 
a random vector with the steady-state distribution of the associated RBM. 
Kaspi and Mandelbaum [23] proved a similar convergence for a closed GJN. 
(This case is quite distinct from the open network setting since the state 
space in the former is compact and, hence, establishing tightness of the 
sequence of re-scaled queue lengths is straightforward.) A very interesting 
result in the context of large deviations theory was established by Majew- 
ski [31]. He showed that, in a feedforward-type multiclass queueing network 
with deterministic service times, the large deviations rates corresponding to 
the network converge to the large deviations rates of the associated RBM 
when the network is in heavy traffic. While these results are suggestive of the 
fact that the stationary distribution of RBM may indeed provide a rigorous 
approximation to that of the underlying GJN, the recent book of Chen and 
Yao [11] mentions that this remains an open problem. Resolving this ques- 
tion would constitute an important stepping stone on route to proving the 
validity of RBM steady-state approximations in other classes of queueing 
networks. 
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Fig. 1. The interchange- of -limits diagram establishing the validity of heavy-traffic 
steady-state approximations. 

The main goal of the present paper is to establish this conjectured cor- 
respondence between the steady-state of GJN and its Brownian approxima- 
tion. In particular, we prove that the stationary distribution of the vector of 
queue lengths in the GJN, under diffusion scaling, converges to the station- 
ary distribution of the associated RBM as the traffic intensity approaches 
unity (see Theorem 8). This result establishes the validity of an "interchange 
of limits" argument. Specifically, by first taking heavy-traffic limits, one ar- 
rives at an RBM approximation from which, by letting time go to infinity, 
one obtains the steady-state approximation to the GJN. What we show is 
that this order of limits can be interchanged. In particular, we prove that the 
same stationary distribution is obtained by first passing to steady-state in 
the GJN, and subsequently applying diffusion scaling and letting the system 
utilization approach one. This argument is represented symbolically in the 
diagram given in Figure 1, where Q n (t) represents the queue length process 
in heavy-traffic, and Z(t) represents the associated RBM. Using the above, 
we also prove that the normalized moments with respect to the stationary 
distribution of the GJN converge to the respective steady-state moments of 
the RBM (see Corollary 2). Finally, we establish the above interchange ar- 
gument for the sojourn times (see Theorem 10), thus proving a steady-state 
version of Reiman's [33] "snapshot principle." 

The main technical subproblem that we solve in order to establish our 
main result is that of tightness of the normalized queue length vector in 
steady-state. Roughly speaking, we show that, as the traffic intensity pj — > 
1 in every station j = 1, ...,«/, the family of scaled queue length vectors 
((1 — Pj)Qj(t))j is tight (see Theorem 7). This establishes the existence of a 
probability measure which is a weak limit of this family. It is then a matter 
of simple argument to show that this measure must be the unique station- 
ary distribution of the associated RBM. Proving tightness of the family of 
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normalized steady-state queue lengths hinges on the following ingredients: 
(i) existence of a pathwise solution to the Skorohod problem (the oblique re- 
flection mapping); (ii) existence of a "fluid-scale" linear Lyapunov function 
for the workload; and (hi) probability bounds for strong approximations of 
the primitive processes, that is, arrivals, services and routing. 

On route to our main result, we derive exponential bounds on the sta- 
tionary distribution of general Markov processes admitting a suitable Lya- 
punov function (see Theorems 5 and 6). While such bounds are excessive 
for purposes of proving tightness of the stationary distribution under dif- 
fusion scaling, they are potentially quite useful for performance analysis of 
GJN. (For recent work that derives tail bounds in the heavy-tailed context 
for GJN, and more general classes of monotone separable networks, see [1].) 
These bounds and their derivation are in the spirit of the results obtained 
by Meyn and Tweedie in several papers and summarized in their book [32]. 
Unlike the aforementioned results of Meyn and Tweedie [32] (see also the 
queueing application in [14]), our work does not use an infinitesimal drift 
characterization of the Markov process in question. Rather, we derive our 
bounds on the basis of the behavior of the Markov process on "large" time 
scales. This approach is more intuitive, supports a more direct application 
of the theory to the GJN context, and circumvents several technical issues. 
In terms of strong approximations, we rely on similar ideas to those used 
by Horvath [21], however, our main focus is on establishing uniform integra- 
bility of diffusion-scaled processes. The blend of methods we employ might 
be of independent interest and, in a parallel ongoing work, we use these 
techniques to derive explicit upper and lower bounds on steady-state queue 
lengths in GJN under both heavy-tailed and light-tailed assumptions on the 
primitives. 

The main technical assumption we impose on the primitives of the GJN 
is the following. Both the distribution of interarrival and service times are 
assumed to have a finite conditional exponential moment generating func- 
tion (a precise definition is given in what follows). We note that "standard" 
heavy-traffic limit theory, including Reiman's theorem, typically only as- 
sumes the existence of a second moment (more precisely 2 + 5 moments) 
for the interarrival and service time distributions. In that sense our condi- 
tions are clearly much stronger than necessary. Nevertheless, these condi- 
tions facilitate the mathematical derivations and, moreover, lead to novel 
exponential tail bounds on the steady-state of the GJN (see Theorem 6). 
Extending our results to the more general case assuming only existence of 
the first two moments is an important direction for future work. Note that 
the RBM corresponding to the heavy-traffic limit of the underlying GJN 
has a steady-state with exponential tails; see [19]. This suggests that the 
existence of exponential moments for the GJN primitives should result in 
a more accurate RBM-based steady-state approximation. To this end, it 
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is worth noting that, even in the simple case of the GI/GI/1 queue with 
heavy-tailed service times, the RBM exponential steady-state is known to 
provide a poor approximation of the Pareto-like stationary distribution of 
the original queue length process. In Section 5 we discuss briefly conditions 
under which we believe our analysis should extend to the case of GJN with 
primitives having Pareto-like distributions. 

The rest of the paper is organized as follows. Section 2 provides back- 
ground and preliminaries on GJN, the Skorohod mapping and heavy-traffic 
approximations. Section 3 derives bounds on the stationary distribution of 
general state-space Markov chains using Lyapunov arguments. Section 4 de- 
tails the main results. Specifically, the main theorem is stated in Section 4.2. 
Concluding remarks are given in Section 5, and all proofs are collected in Ap- 
pendix. 

2. The queueing network model and its heavy-traffic approximation. The 

set up and notation follow closely that in [11], Chapter 7. We use / to denote 
a J x J identity matrix, and e the J-dimensional vector of ones. (All the 
vectors are assumed to be column vectors, unless otherwise stated.) Put e l 
to be the ith unit vector in W J . The transposed of a matrix P (vector e) 
will be denoted P' (e'). For every z € M. J , the norm ||z|| corresponds to the 
L\ metric: ||z|| = X)i<j<j \ z j\- Given a random variable X and a probability 
measure z/, X ~ v means X is distributed according to v. For a sequence of 
probability measures v n , n = 1, 2, . . . , v n =>■ v is used to denote weak conver- 
gence of v n to a limit probability measure v. A similar notation is used to 
denote weak convergence of a sequence of random variables or processes. 

2.1. Generalized Jackson networks: model description and probabilistic 
assumptions. 

Description of the network and probabilistic assumptions. The network 
consists of J stations denoted for simplicity by 1, 2, . . . , J. Jobs arrive from 
external sources to each station according to independent renewal processes 
with interarrival times given by dj(0), Oj(l), «j(2), . • • for j = I,..., J. The 
interarrival times aj(l), a, (2), . . . are assumed to have a common distribution 
Faj for j = 1, . . . , J. The first interarrival time dj(0) is assumed to have the 
same distribution function, only conditioned on the event cij(O) > clqj, where 
a o,j > measures the time elapsed since the last arrival prior to time t = 0. 
This is a given deterministic or random value and is considered to be part of 
the data. Thus, the distribution of the arrival process is completely specified 
by a pair (Faj^oj)- Put ao := (ao,i> ■ • ■ > «o,j)- Also, let ctj = 1/E[a_j(l)] 
be the external arrival rate into station j, and put a = (a%, . . . ,ctj) to 
be the vector of external arrival rates. It is not excluded that ctj = for 
some stations j, meaning there is no external arrival into these stations. 
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We denote by J = {j :aj > 0}. Let Aj(t) = sup{n:J2o<i< n a j(l) — denote 
the corresponding counting renewal process. For every ti > t\ > 0, we put 
Aj(ti,t2) ■= Ajfa) — Aj{t\). The coefficient of variation var[a|(l)]a| is de- 
noted by (? a p for j e J . 

Denote by Vj(0),Vj(l), . . . the sequence of service times received by jobs 
in station j and Vj(n) = J2o<i< n v j(0- The service times Vj(l),Vj(2), . . . , are 
assumed to be i.i.d. with a distribution function Fsj- As in the case of in- 
terarrival times, the first service time Vj(0) is assumed to follow the same 
distribution, only conditioned on Vj(0) > voj. Here voj > is a (possibly 
random) value which stands for the cumulative service time elapsed for the 
job currently in service in station j at time t = 0. This value is considered a 
part of the data. For consistency, we assume that Vqj > only if the queue 
length in station j at time is positive. We let Sj(t) = sup{n : Vj(n) < t} de- 
note the associated renewal process and let rrij = M[vj(l)] denote the mean 
service time in station j. Then, fij = 1/mj is the service rate in station j, 
and fj, = (m, . . . ,fj,j) is the vector of service rates. We let M be the diag- 
onal matrix with m = (mi, . . . ,mj) as diagonal entries. The coefficient of 
variation var[u|(l)]^| is denoted by (? s •, for j = 1, . . . , J. 

We assume that both the interarrival and service times have a finite mo- 
ment generating function in some neighborhood of the origin. Moreover, to 
handle the residual service and interarrival times, we assume the following 
condition holds: there exists 9* > such that, for every j, 

(1) sup E[exp(0*(aj(l) - z))K(l) > z] < oo, 

(2) sup E[exp(0*(uj(l) - z))\vj(l) > 

2GM+ 

where j G J in (1) and j = 1, 2, . . . , J in (2). In particular, E[exp(0*aj(l))] 
and K[exp(9*Vj(l))} are finite as well. 

The routing decisions at each station are assumed to be of Bernoulli type 
and parameterized by a sub-stochastic J x J matrix P= (pjk)i<j,k<j- The 
entries of P satisfy pij > 0, for every i,J2jPij — 1 an d the spectral radius of 
P is strictly less than unity. For each station j, let i/j 3 = {ip J (l),^ (2), • • •} be 
an i.i.d. sequence of routing decisions with common distribution Pj, where 
Pj is the jth column of P. Specifically, = e k ) =Pjk an d W(ipi(l) = 

0) = 1 — J2i<k<jPjk, the latter indicates the fact that the completed job 
leaves the network. Let i? J (0) = and B?(n) = Ekkh^O)- Then the fcth 
component, i?i(n), of i?- ? (n) is the total number of jobs (out of n) which, 
upon completing their service at station j, are routed to station k. 

Let A = (Ai, . . . , Aj) denote the unique solution to the traffic equation 
A = a + P'A. That is, A = [I - P'^a. Then p = M~ 1 A = M' 1 ^ - P']a = 
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(Pj)i<j<j is the vector of traffic intensities. We let p* = m&Xj pj denote the 
bottleneck traffic intensity. We will primarily consider the case when 

(3) p* < 1. 

This is also referred to as the stability condition for GJN, a terminology 
which is justified by Theorem 2 below. 



Systems dynamics. The number of jobs which are either in service or 
waiting in station j at time t is denoted by Qj(t); adhering with standard 
terminology, we refer to Qj = (Qj (t):t> 0) as the queue length process at 
station j. Put Q = (Qi, ... , Qj). Let Bj(t) denote the cumulative time that 
station j was busy processing work during the interval [0, t) . In particular, 
< Bj(t) < t,B(0) = 0. Put Bj = (Bj(t):t > 0) and let B = (B u . . . , Bj); 
we refer to this as the busy time process. Let Ij(t) = t — Bj(t) denote the 
cumulative idle time at station j in the interval [0, t]. Function Ij = (Ij(t) : t > 
0) describes the idle time process in station j. At every station j the jobs 
are assumed to be processed in First-In-First-Out (FIFO) fashion. 

The following equations that describe the system dynamics tie together 
the processes introduced above: 

(4) Q j (t) = Q j (0) + A j (t)+ £ R^SkiBkm-S^it)) 

l<k<J 

and 

(5) Ij(t)= [\{Q j (s) = 0}ds, 

Jo 

for all t > 0, j = 1, 2, . . . , J. We consider the centered process X = (Xj (t):t> 
0) with 

Xj (t) = Qj (0) + ( aj + ViPij ~ H ) * + ( A J (*) - 

V l<i<J / 

(6) + E PijiSiiBM-^B^-iSjiBjit^-fijBjit)) 

l<i<J 

+ j2 (Riismm-pijSiiBiit))), 

l<i<J 

j = 1, 2, . . . , J, and the process Yj = (Yj(t) :t > 0), where 

(7) Yj{t)= N Ij{t)= N {t-Bj{t)). 

The process Xj is often referred to as the "free process" which summarizes 
the net input of jobs to station j. The process Yj summarizes system idleness 
in normalized units. Alternatively, this process "regulates" the queue lengths 



STEADY-STATE APPROXIMATIONS IN OPEN QUEUEING NETWORKS 9 



from taking negative values by increasing whenever the queue length process 
is zero. The dynamics in (4)-(5) can now be represented in the following 
succinct form: 

(8) Q(t) = X(t) + [I-P']Y(t)>0, 

(9) y(0) = an d y(t) is nondecreasing in t, 

(10) / Qj (r) dYj (r) = for all t > 0, 1 < j < J. 
Jo 

The equations (8)-(10) describe the system dynamics (Q,Y) as the solu- 
tion to the Skorohod reflection problem associated with the GJN. A formal 
definition of the Skorohod problem is given in Section 2.2 below. 

A Markovian representation. We now provide a Markov process repre- 
sentation for the underlying queue length processes. The queue length pro- 
cess Q is clearly not Markovian because of the residual interarrival and ser- 
vice times. To create a Markovian structure, we augment the state descriptor 
as follows. For each time t > and j £ J ', we let dj (t) denote the time elapsed 
since the most recent arrival to station j that occurred prior to time t. That 
is, dj{t) := t — Y^n<Aj{t) a j{ n )-> an d we put a(t) = (fii(i), . . . ,dj(t)), where, 
by definition, dj(t) = for j ^ J for all t. Similarly, Vj(t) denotes the time 
elapsed since the commencement of the most recent service prior to time t 
at station j, or if the server is idle at time t, this variable is equal to zero. 
Put v(t) = (vi(t), . . . ,vj(t)). The extended process Q = (Q(t) : t > 0) with 
Q(t) = (Q(t),a(t),v(t)) is a Markov process. That is, for every t\ > t > 
and every (z,a,v) G x M^ 7 , the distribution of Q{t\) conditioned on 
Q(t) = (z,a,v) is independent of Q(t') for all t' < t. (In fact, the process 
Q is known to be a strong Markov process; see [11] for further details.) By 
definition, Q(0) = (Q(0),a(0),#(0)) = (Q(0),ao,i>o)- For brevity, we denote 
the extended state space of the process Q by X := Zi x R 2 / . The probability 
distribution of the strong Markov process Q(t) for every t > is completely 
specified by (Q(0), do, i>o, Faj, Fsj, P)- We denote this 6-tuple by H and call 
it the parameter vector for the corresponding generalized Jackson network. 
For brevity, we often refer to the underlying GJN simply as H. A probabil- 
ity distribution ir on S is defined to be stationary if Q(0) ~ ir implies that 
Q(t) ~ 7r for all t > 0. [The stationary distribution for the Markov process 
(Q(t) :t>0) is then constructed in the usual manner from n and the tran- 
sition kernel.] When Q(t) ~ tt we say that the GJN is in steady state with 
stationary distribution having ir as its marginal; this distribution is not 
necessarily unique. For future purposes, we simply use ir when making ref- 
erence to the stationary distribution of the Markov process in question. We 
will also say that the network is stable if there exists at least one stationary 
distribution ir. 
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2.2. The Skorohod mapping, fluid model and stability of GJN. This sub- 
section reviews some results concerning the Skorohod problem and fluid 
analysis of GJN that will be useful in what follows. The presentation is 
based on Chapter 7 in [11]. 

Skorohod problem and the oblique reflection mapping. Let £>[0,oo) de- 
note the space of right-continuous functions with left limits taking values in 
~R J . Since all limit processes considered in this paper have continuous sample 
paths, it is not necessary to introduce the Skorohod topology, and it suffices 
to consider the space T> endowed with the uniform topology; see [6] for fur- 
ther discussion. We now cite a general result, originally derived in [18] in 
the context of reflected Brownian motion, that supports the representation 
of GJN dynamics given in (8)-(10); the version below is adapted from [11], 
Theorem 7.2. 

Theorem 1 (Oblique reflection mapping). For every function x £T>[0, oo) 
such that x(0) G R+, there exists a unique pair of nonnegative functions 
(y, q) = (^i(x), ^ 2 (x)) £ T>[0, oo) such that yj(t) is a nondecreasing function 
for every j = 1, 2, . . . , J, yj(0) = 0, and (x, y, q) jointly satisfy 

(11) q(t)=x(t) + [I-P']y(t), 

(12) / q j {T)dy j {T) = Q for all t> and j = 1,2, ... ,J. 
Jo 

Moreover, the mappings *$>i, *$> 2 '■ f[0, oo) — > T>[0, oo) are Lipschitz continu- 
ous, in particular, for every Xj,^ € 22 [0, oo) and every t>0, 

(13) sup ||<7i (r) -q-i{r) \\ <R\ sup \\xx{t) - x 2 (t)\\ 

0<T<t 0<T<t 

and 

(14) sup \\yi(T) - y 2 (r)\\ < R 2 sup ||xi(r) - x 2 (r)||, 

0<T<t 0<T<t 

where R\,R 2 > depend only on the matrix P. 

Remark 1. The processes q and y are referred to as the reflected pro- 
cess and the regulator process, respectively. For further details and explicit 
identification of the Lipschitz constants R\,R 2 , see, for example, the recent 
book by Whitt [39], Chapter 14.2.3. 

Remark 2. The existence and uniqueness of the pathwise solution to 
the Skorohod problem hinges on the matrix P being substochastic. This, in 
turn, implies that the reflection matrix I — P' is Minkowskii, that is, has all 
positive diagonal elements, nonpositive off-diagonal elements, and inverse in 
which all entries are nonnegative. 
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The fluid model and stochastic stability. Given any z £ M+, let x z (t) = 
z + (a - (I - P')p)t and let (y, q) = ^(z))- Tne triplet (x,y,q) is 

called a /Zmd model of the GJN. Let w = e'[I — P'] -1 . Note that each coor- 
dinate 101, . . . , w j, of the vector 

w = e'[I-P'}- 1 

= e'[I + P' + (P') 2 + (P'f + ---] 
is greater than one. Consider the process 
(15) (w'q(t):t>0). 

It is well known that w'q(t) is Lipschitz continuous (see [11]). The process 
(w'q(t) :t>0) can be interpreted, essentially, as the total amount of "fluid 
work" present in the system at time t. An important property of the fluid 
model is the amount of time needed to drain existing work from the system 
when all exogenous inflows are "turned off" at time t. For the fluid model 
of a GJN network, we have the following result, which can be found, for 
example, in [11]. 

Proposition 1 (Fluid stability of GJN). The vector valued function 
q(t) is differentiable everywhere on R+ except for a set of points which has a 
Lebeasgue measure zero. Moreover, w'qit) < — mini<j<j/ij(l — pj) whenever 
q(t) 7^ and q(t) is differentiable. As a result, if condition (3) holds, then 
for any initial fluid level q(0) = z, there exists 



mmi<j<j fij(l - pj) 
such that q(t) = for all t>r. 

The function w'q(t) serves as a Lyapunov function for the GJN network. 
In particular, the result articulated above asserts that the fluid model is sta- 
ble, as it drains in finite time from any initial fluid level. Given the results of 
Dai [12] and Stolyar [37] linking stochastic and fluid stability, Proposition 1 
can be used to establish stability of the underlying GJN. To this end, Sig- 
man [35], Down and Meyn [15], Dai [12] and Dai and Meyn [14] prove the 
following theorem. 

Theorem 2 (Stochastic stability of GJN). Suppose condition (3) holds 
for the GJN S. Then, there exists a stationary distribution ir for the Marko- 
vian process Q{t). Under additional technical regularity (e.g., Assumptions 
A1-A3 in [14],), the stationary distribution is unique and for every starting 
state 0(0) = (z,a,v) G X, the distribution of Q(t) converges to tt as t — > oo. 

The Lyapunov function w'q(-) and the fluid model of the GJN will play 
an important role in deriving our main results in what follows. 
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2.3. Heavy-traffic limits and reflected Brownian motion. 

The GJN in heavy-traffic. We now introduce the generalized Jackson 
network in the parameter regime corresponding to heavy-traffic assumptions. 
We consider a sequence of networks with traffic intensities pj,n = 1,2, ... , 
each approaching unity at the same rate as n — > oo. To arrange for that, we 
first consider a network such that the traffic intensity pj = 1 in every station 
j, and denote this network by H. (At this point the reader should recall that 
H encodes the 6-tuple given in Section 2.1 that parameterizes the GJN.) In 
particular, X = [I — P']~ 1 a = p. Fix a vector of strictly positive constants 
kP = • • • Then, we consider a sequence of networks {H n } in which: 
(i) the service times Vj (as well as the residual service times) and routing 
decisions ip 3 have the same distribution as in the original network S; (ii) 
the interarrival times for the nth network in the sequence are defined to be 
a"j = a,j(l — Kj/y/ri), and the residual interarrival times remains unchanged. 
Hence, the vector of arrival rates is given by a n = (a™, . . . , a"). As a result, 
the vector of effective arrival rates is X n = [I — P / ]~ 1 a n , which we write as 
[/ - P']- l a - n- l l 2 [I - P'}- l K a = A - n- l ' 2 K, where «:=[/ — P']- 1 ^ = 
(ki,..., kj). Then 

(17) fi = (l-Kj/yfii) for 3 = 1,..., J. 

Equivalently, >Jn(\ — p™) = Kj for j = 1, . . . , J. We also obtain for the bot- 
tleneck traffic intensity p* n = (1 — k* j \fn), where k* = minKj. The GJN in 
heavy-traffic refers to the sequence of generalized Jackson networks {H n }, 
and the corresponding sequences of processes are denoted A n (t), B n (t), 
I n (t), X n (t), Q n (t), Y n (t). Note that the processes S(t) and the routing 
processes do not depend on n, so we use the notation S(t),R J (l) for these 
elements in the network H n . The sequence of stationary distributions cor- 
responding to H n will be denoted by {ir n } henceforth. These distributions 
exist by Theorem 2 due to the fact that condition (3), p* n = m&Xj p™ < 1, 
holds for all n. When the stationary distributions are not unique, we denote 
by {7r n } an arbitrary sequence of stationary distributions corresponding to 

Reflected Brownian motion in the orthant. The importance of the heavy- 
traffic regime stems from the ability to approximate the underlying queue 
lengths and workload processes by a reflected Brownian motion (RBM) 
which we define below. 

Definition 1 (Reflected Brownian motion in the orthant). Let z £ M J , 
P be a substochastic J x J matrix, and let W = (W(t) :t > 0) be a J- 
dimensional Brownian motion with initial state z, drift vector j3 S M" 7 and 
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covariance matrix T 6 x . Then, the reflected Brownian motion (RBM) 
with parameters {(5, T,I — P') and initial state z is the unique solution to the 
Skorohod problem (11), (12) with input process W. In particular, the RBM 
Z = (Z(t) :t > 0) is given by Z = ^i(W), and the nondecreasing regulator 
process Y = (Y(t):t > 0) is given by Y = * 2 (W). Thus, the pair (Z,Y) 
satisfies 

(18) Z(t) = W(t) + [I-P']Y(t), 

(19) /* Zj (s) dYj (s) = for all t > and j = 1, 2, . . . , J. 
Jo 

The RBM process Z has the following behavior: in the interior of the non- 
negative orthant M.+ evolves like a Brownian motion with drift (3 and diffu- 
sion matrix T; and when the process hits the jth face, Fj = {i£ : Xj = 0}, 
it is reflected instantaneously toward the interior of the orthant in the di- 
rection determined by the jth column of the matrix I — P' . Roughly speak- 
ing, the jth component of the regulator process Y measures the cumulative 
"effort" in maintaining the Brownian motion confined to the nonnegative 
orthant. 

The stationary distribution of the RBM (when it exists) is one of the 
primary focuses of the present paper. To that end, the following result, 
established in [19], will be important in what follows. 

Theorem 3 (RBM stationary distribution). The RBM (j3,T,I - P') 
possesses a stationary distribution tt^bm if and only if [I — P']^ 1 j3 < 0. When 
this distribution exists, it is unique. Moreover, for every initial distribution, 
the marginal distribution of RBM converges weakly to 7Trbm as time diverges 
to infinity. 

Heavy-traffic scaling and diffusion limits. For the sequence of GJN, {H n }, 
defined earlier, put 

X n (t)=X n (nt)/V^, 

(20) 

Q n (t) = Q n (nt)/VE. 

The above correspond to a diffusion scaling of the free process X and the 
queue length process Q. In particular, the normalization in (20) corresponds 
to accelerating time and compressing space in accordance with the functional 
central limit theorem (FCLT). The following theorem, due to Reiman [33], 
is one of the most fundamental results in the theory of diffusion approxima- 
tions of queueing networks. The theorem was originally proven under the 
assumption that Q n (0) = 0, however, the more general result stated below 
also holds true (see, e.g., [11]). 
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Theorem 4 (Process limit weak convergence of GJN to RBM). Con- 
sider a sequence of GJN, {H n }, with p n satisfying (17), and assume that 
Q n (0)/T,/n Z(0), where Z(0) is distributed according to some probabil- 
ity measure ttq on Ri. Then, for every t > 0, the sequence of processes 
(Q n (nt')/y/n:0 < it < t) converges weakly to the RBM (Z(t'):0 < it < t) 
with respect to the topology of uniform convergence in T>[0, t] . The RBM has 
initial distribution Z(0) ~7To, and it has parameters ((3,T,I — P'), where the 
drift vector is given by (3 = —(I — P')M~ l k, and the entries of the covariance 
matrix T are given by 
J 

Tk£ = ^2vj[Pjk(fiki-Pj£) + c %j(Pjk - 5jk)(Pj£ - Sji)] +a k c 2 a k 5 M , 

where 5jk = 1 if j = k and is zero otherwise. 

Remark 3. The limiting covariance structure can be inferred from (48) 
in Section 5; for further details, see [11], Theorem 7.29. 

Remark 4. Note that the expression for the RBM drift vector implies 
that [I — P'\~ 1 j3 = —M~ 1 k < 0, which is the condition for the existence and 
uniqueness of the stationary distribution of RBM, and corresponds to the 
condition p* n < 1 in H n . 

Remark 5. An important point concerning the RBM approximation is 
that it can give rise to closed form solutions for its stationary distribution, 
while the underlying GJN may be intractable in that regard. Harrison and 
Williams [19, 20] provide necessary and sufficient conditions under which the 
stationary distribution of RBM is separable with exponentially distributed 
marginals. An example that illustrates this is given in Section 4. 

3. Lyapunov functions and steady-state bounds for Markov processes. 

3.1. Background and definitions. Let S = (H(£) - t > 0) be a continuous 
time Markov process defined on a complete metrizable state space X, which, 
for our purposes, is Euclidean space, equipped with Borel <r-algebra B. For 
the special case of generalized Jackson network, H(t) = Q(t) = (Q(t),a(t),v(t)) 

Let Pa; be the probability distribution under which P(3(0) = x) = 1, and 
put K x [-] :=E[-|H(0) = x], that is, the expectation operator w.r.t. the prob- 
ability distribution ¥ x . Let ¥ n denote the probability distribution under 
which S(0) is distributed according to n, and put E„-[-] to be the expecta- 
tion operator w.r.t. this distribution. A probability distribution ir defined 
on B is said to be a stationary distribution if for every bounded continuous 
function / : X — > R, 



(21) 



E T [/(S(0))]=E 7r [/(H(t))], 
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for all t > 0. The following definition plays an important role in our analysis. 

Definition 2 (Lyapunov function). A function : X — > M + is said to be 
a Lyapunov function with drift size parameter —7 < 0, drift time parameter 
to > and exception parameter K, if 

(22) sup {E**(3(t )) - *(x)} < -7. 

xe.#:*(a;)>.K 

A function $ : A* — ► M + is defined to be a geometric Lyapunov function with 
a geometric drift size < 7 < 1, drift time to an d exception parameter K, if 

(23) sup {($(x))- 1 E a; $(H(to))}<7- 

xeX : <Z>(x)>K 

Given any function <3? : X — > K+ and t > 0, let 

(24) <P(t) = snp^ 1 (x)E x ^(E(t)). 

It is not excluded that (j>(t) = 00, however, the interesting case of course is 
when <p(t) is finite. Given 6 > 0, let 

(25) Li(0,t) = supE x [exp(^($(S(t)) - $(x)))] 
and 

(26) L 2 (0,t) = supE x [($(H(t)) - $(x)) 2 exp(0($(H(t)) - $(x)) + )] 
for t > 0. 

3.2. Bounds on the stationary distribution. Lyapunov functions provide 
a useful tool for obtaining bounds on stationary distributions of Markov 
chains and processes. In this subsection we are set to do just that. The 
bounds we derive are very similar to the ones obtained in various forms and 
under a variety of assumptions by Meyn and Tweedie, summarized in their 
book [32]. 

Theorem 5. Suppose the Markov process (H(t) :t > 0) possesses a sta- 
tionary distribution ir and suppose $ is a geometric Lyapunov function with 
parameters 7, to, K . Then, 

4>(t )K 



(27) ^[$(5(0))] < 



1-7 



Note that, in the theorem above, it is not assumed that the stationary 
distribution tt is unique. The claimed bounds hold for every stationary dis- 
tribution 7T. The bound of course is only meaningful when (f>(to) is finite. 

We now apply Theorem 5 to establish bounds on the tail of the stationary 
distribution. 
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Theorem 6. Let (B,(t) :t>0) be a Markov process with a stationary dis- 
tribution tt, and suppose <3? is a Lyapunov function with parameters -y,to,K. 
Assume, in addition, that there exists 9 > such that 

(28) 0L 2 (0,i o )<7- 

Then exp(0<I>(-)) is a geometric Lyapunov function with geometric drift size 
parameter 1 — 70/2, drift time parameter to and exception parameter e eK . 
Consequently, for every s > K, 

(29) P*(*(E(0)) > s) < (1 - 70/2)-% (0, t ) exp(-0(s - K)). 

Remark 6. At first glance, the division by 1 — 70/2 above may seem 
erroneous since, in principle, this expression can be negative. We now show 
that condition (28) rules this out. Prom the definition of L 2 (v), we have 
that 

L 2 (0,i o ) > su P E x [($(H(t ))-$(x)) 2 ] 
> sup (E^ [<I>(£(to)) - $(x)]) 2 




where (a) follows from Jensen's inequality and (b) follows from the definition 
of 7, in particular, 

sup E a .[*(H(t )) - < -7- 

x:&(x)>K 

The condition (28) then implies that 70 < 1. 



Remark 7. The upper bound in (29) can be optimized by selecting 
> to be the largest possible value satisfying (28). Of course, if the largest 
such is equal to zero, then the bound is trivial. 

4. Main results. Throughout this section we consider a sequence of gen- 
eralized Jackson networks H n in heavy traffic, and the associated Markov 
process Q n (t) = (Q n (t) , a n (t) , v n (t)) introduced in Section 2; the latter two 
coordinates denote the vectors of residual interarrival times and service 
times, respectively. 

4.1. Tightness of stationary distributions. In this section we establish 
tightness of the sequence of stationary distributions corresponding to the 
normalized state process {Q n (nt) j \/n} . Key to this result are the tail bounds 
derived for general Markov processes in Section 3. To appeal to these results, 
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we consider the workload function denned in (15), whose key properties are 
given in Proposition 1. The next result establishes that this workload func- 
tion has the requisite attributes to serve as a Lyapunov function for the 
sequence of GJN's in heavy-traffic. 

Proposition 2. There exist constants toi c o > which depend only on 
the parameters of the network E and are independent of n, such that, for all 
sufficiently large, n > 1, 

(30) sup {E[w'Q n (nto)\Q n (0) = (z,a,v)]-w'z}<-y/n. 

(z,a,v) : w' z>co^/n 

In addition, there exists 0q > which depends only on S such that 

limsup sup E[exp(n- 1/2 {w'Q n (nt )-w'z) + )\Q n (0) = (z,a,v)] 

n^oo ( Zt a,v)eX 

(31) 

< oo 

and 

limsup sup n^ 1 K[(w' Q n (nto) — w' z) 2 

n-^co (z,a,v)eX 

(32) x exp(n- 1/2 e (w'Q n (nt ) - w'z) + )\Q n {0) = (z,a,v)} 
< oo. 

The immediate implication of the above result is the following. 

Proposition 3. For all sufficiently large n, the function $>(z, a, v) = w'z 
is a Lyapunov function with drift size parameter —y/n, drift time parameter 
nto and exception parameter coy/n. In addition, 

(33) limsup Li(9q / y/n, nto) < oo, 

n-^oo 

(34) limsup— L2(0q/ yfn,nto) < oo. 

n— >oo n 

Armed with this result, we now use Theorem 6 to establish the sought 
tightness of the sequence of stationary distributions corresponding to the 
sequence of normalized GJN's. 

Theorem 7. There exist constants C\,c\ > which depend only on S 
such that, for, all sufficiently large n, the sequence of stationary distributions 
7r n of the networks {B n } satisfies 

(35) P 7r »(n -1/2 u; / Q ri (0) > s) < dexp(-c lS ), 
for all s > 0. 
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Since 1 — p* n = k* /y/n, the theorem implies that each individual queue 
Qj'it) behaves in steady state as 0( 1 _ 1 pm ) as p* n — > 1. 

Recall that a sequence of real valued random variables X n ,n G N, is said 
to be tight if, for every e > 0, there exists sufficiently large constant C > 
such that sup n P(|X n | > C) < e. Refer to [6] for more detailed definitions 
and applications of this important concept. The following corollary is a con- 
sequence of the theorem above. 

Corollary 1. Let Q n (0) = (Q n (0),a n (0),?) n (0)) be distributed accord- 
ing to 7r n . Then the sequence of I'' -valued random vectors {Q n (0)/y/n} is 
tight. 

4.2. Interchange of limits. We are now in position to prove the sought in- 
terchange of limits indicated in the schematic diagram in Figure 1 . Let 7r n de- 
note the stationary distribution of the rescaled process Q n (t) = Q n (nt)/^/n, 
where, as above, (Q n (t) :£ > 0) denotes the queue length process for the GJN 
with traffic intensity p n . With respect to 7r n , we have, by Theorem 2, that, 
under some mild assumptions for each fixed n, the distribution of Q n (t) 
converges weakly to tt u as t — > oo. On the other hand, Reiman's diffusion 
approximation stated in Theorem 4 asserts that starting with the origi- 
nal GJN and taking heavy-traffic limits, that is, letting p n f 1 as n — > oo, 
the normalized queue length process (Q n (t) :t>0) converges weakly to an 
(J3, r, J - P)-RBM (Z(t) : t > 0) identified in Section 2. Subsequently, letting 
t — > oo, one obtains a steady-state approximation for the original network 
given by the distribution vtrbm- We now prove that the above limits can 
be interchanged, that is, first taking the limit as t — > oo and letting the 
GJN reach steady-state, then taking heavy-traffic limits; the resulting pro- 
cess converges to a weak limit and this limit is 7Trbm- This validates the 
steady-state heavy-traffic approximation which is built on the distribution 
of the approximating RBM. 

To prove the interchange of limits recall the assertion of Corollary 1. 
In the space of probability measures on M. J endowed with topology of weak 
convergence, the tightness result of this corollary implies that there exists at 
least one limit point, tt* , for the sequence {tt"}, and this limit is a probability 
measure. Our main result stated below asserts that such a limit is unique 
and is identical to the unique stationary distribution of the associated RBM. 

Theorem 8. The sequence {tt 71 } converges weakly to 7Trbm> where ttrbm 
is the unique stationary distribution of the ((3,T,I — P)-RBM. 

Thus, we have that the probability distribution of {Q n (0)/\/n} converges 
weakly to 7Trbm when Q n (0) is distributed according to tt 71 . This establishes 
the validity of the commuting diagram depicted in Figure 1. 
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Given our assumptions on the primitives, in particular, the existence of 
a (conditional) moment generating function in the neighborhood of zero, it 
stands to reason that the queue length process under diffusion scaling should 
be uniformly integrable. Moreover, we expect the interchange argument to 
hold for all moments of the normalized queue length process. The following 
result asserts that this is indeed the case. 

Corollary 2. There exists a vector 6 E with all coordinates being 
strictly positive, such that the family of random variables {exp(^'Q n (0)/y / n)} 
is uniformly integrable, where Q n (0) ~ 7r™ for each n= 1,2, Conse- 
quently, for all p> 1, 

[n-^ 2 Q n (0)y E nRBM [Z(0)] p as n - oo, 

where vtrbm is the stationary distribution of the RBM (Z(t):t > 0) with 
parameters ((3,T, I — P'). 

We note that the above result is particularly useful given that current nu- 
merical algorithms for analysis of the steady-state of RBM, see [13] and [36], 
are mostly effective in the computation of stationary moments (as opposed 
to the full distribution). 

4.3. Extension to sojourn times. Given a job which arrives externally 
into station j E J ', we will say that this job has the visit vector h E if it 
visits station i precisely hi times before leaving the network. Of course, not 
every visit vector is feasible. That is, there may not exist a job with visit 
vector h for all h E Zi. For example, in a tandem queueing system with two 
stations and no feedback, only h = (1, 1) is feasible. As in [33], we let Dj^t) 
be the sojourn time of the first job which arrives into station j after time 
t and has a visit vector h. Given a state Q(t) = (z,a,v) E X of the Markov 
process at time t, the distribution of Da ^(t) is independent from Q(t'),t' < t. 
Thus, we may consider the random variable Dj t h(z,a,v) which is a function 
of the state (z, a, v). The following result is then an immediate consequence 
of the above observation. 

Proposition 4. Let n be the stationary distribution of the Markov pro- 
cess Q(t). Then for every station j and feasible visit vector h E r Li r , the 
process Dj t h(t) := Dj^{Q(t)) is stationary. 

In the queueing system E n let D™ h (t) denote D™ h (Q n {t)). The weak con- 
vergence of the processes D^(nt) j \fn defined in the queueing network S n to 
an RBM was established by Reiman [33]. The result is stated there under 
the assumption that the queueing network is initially empty Q n (0) = 0. The 
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techniques used by Reiman allow one to obtain a stronger result, which we 
will use in the present paper. We do not give a full proof of this stronger 
version, but just highlight the necessary adjustments to the basic arguments 
given in [33]. 

Theorem 9. Consider a sequence of GJN, {H n }, in heavy traffic. Sup- 
pose Q n (0)/y/n converges weakly to some limiting distribution ttq. Then for 
every station j and every feasible h, the process (D™ h {t') / \fn : < t' < nt) 
converges weakly to (h! MZ(t') :0 < tf < oo), where Z(t) is the (f3,T,I — P')- 
RBM with Z~tt . 

We now state our main result for the sojourn times in steady-state. The 
proof is immediate by combining Theorems 8 and 9. 

Theorem 10. Consider the sequence of GJN, {H n }, in heavy-traffic. 
Suppose Q n (0) ~ vr n . Then for every station j and every feasible h G the 
convergence 



holds, where Z(0) ~ 7Trbm- 

4.4. Product form RBM. As established in [19], in special cases described 
below the stationary distribution of RBM has an explicit product form ex- 
ponential distribution. Our main result, Theorem 8, can then be used to 
validate closed form steady-state approximations for certain classes of GJN 
in heavy traffic. In this subsection we illustrate this application of our main 
theorem. 

The following result follows from a more general theorem on RBM in 
steady state, established by Harrison and Williams [19]; see also [11], The- 
orem 7.12. 

Theorem 11 (Product form RBM stationary distribution). Given a 
((3,T,I - P')-RBM, suppose [/ - P']" 1 /? < and, in addition, 



where D is the diagonal matrix with elements (1 —pjj),l <j<J, and A 
is the diagonal matrix with elements Tjj, 1 < j < J. Then the stationary 
distribution vtrbm has a product form density 




^//Gur^o),...,^ 1 ^)) 



(36) 



2T=[I- P')D~ l k + AD^II - P] 



f(zi,...,zj)= Y[ Vjew(-rijZj) 

1<3<J 

where r] = -2A- l D[I - P']" 1 /?. 



ZJ ) g K: 



+' 
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It can be checked that condition (36), which is also referred to as the skew 
symmetry condition, holds when c\ ■ = c\ ■ = 1. In particular, this includes 
the classical Jackson network case. Consider the GJN, {H n }, in heavy traffic, 
and recall that (3 = — [I — P']M _1 k with k = {k\ kj). Using this, we have 
rjj = 2T~- (1 —pjj)/j,jKj. Combining with Theorem 8, we obtain the following 
theorem, which gives an explicit expression for the limiting steady-state 
distribution of GJN in heavy traffic. 

Theorem 12. Consider a sequence of GJN, {S n }, in heavy traffic. Sup- 
pose condition (36) holds. Then the sequence of stationary measures n n sat- 
isfies 

lim ¥ n n I — =- >z\= J] ex pi 2r ; /(l l>jj)lijf'j\i)- 
v i<j<J 

for every z 6 Hi_ . 

Recalling that pj = 1 — Kj/y/n, we can informally express the assertion of 
the above limit theorem as 

ffV((i - Pi)QM > Vi) ~ n ex P (-2ry/(i - p^njzj), 

i<j<j 

when pj ml, for all j = 1, . . . , J, where tt is a stationary distribution of GJN. 
Thus, for certain classes of GJN for which the skew symmetry condition 
holds, we obtain a simple approximation to the steady-state on the basis of 
the product form RBM stationary distribution. 

5. Concluding remarks. We considered a generalized Jackson network 
(GJN) in heavy traffic and proved that the stationary distribution of the 
GJN under diffusion scaling converges to the stationary distribution of the 
associated reflected Brownian motion, thus resolving the open problem of 
so called "interchange of limits" for this class of open queueing networks. 
There are several interesting questions which remain unresolved and are 
worth pursuing in future research. 

1. The bounds on the stationary distribution of the queue lengths in GJN 
obtained in the present paper are not explicit, that is, their dependence on 
the stochastic primitives is not articulated in full detail. The techniques used 
in the paper can be applied, in principle, to obtain more explicit expressions. 
We are currently pursuing such results in ongoing work, deriving explicit 
upper and lower bounds on the steady-state distribution of a GJN. 

2. The assumptions adopted in the present paper are certainly not the 
tightest possible. In particular, it should be possible to carry out a similar 
analysis when one only assumes finite 2 + 5 moments, for any 5 > 0, for 
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the interarrival and service times. (These conditions are very close to being 
necessary and sufficient for the functional central limit theorem to hold, 
and thus for proving heavy-traffic process limits.) The Lyapunov approach 
we develop in the current paper can be adjusted to provide polynomial tail 
bounds under such assumptions; in essence, a pth. moment of the primitives 
would yield existence of the (p— l)st moment of the stationary queue length 
(see, e.g., [14]). This would be used in conjunction with strong approximation 
coupling inequalities that are available under polynomial moment conditions 
on the primitives, the weakest being the existence of a 2 + 5 moment. In 
the present paper we assume exponential moments and consequently obtain 
exponential tail bounds on the steady-state distribution of the GJN. 

3. The main results in our current paper hinge on the use of a simple 
linear Lyapunov function. There is a long history of applying Lyapunov 
function methods for purposes of performance analysis, specifically for mul- 
ticlass queueing networks; see [2, 3, 4, 5, 26, 29, 30]. It should be noted, 
however, that in almost all the existing literature, Markovian-type queue- 
ing networks are considered. The extension of these bounds to general type 
queueing networks has not been explored to date. 

4. This paper presents what is, in our view, a very interesting link be- 
tween diffusion approximations of queueing networks and Lyapunov function 
methods. The latter have been traditionally used to establish stability and 
obtain performance bounds for queueing networks, as well as diffusion limits 
that correspond to the underlying queueing model (see [16] for an example 
of the latter). To the best of our knowledge, such methods have not yet 
been applied to queueing networks to establish properties of their diffusion 
approximations (e.g., to ascertain tightness of the sequence of steady-state 
distributions, as is established in the current paper). The proof technique in- 
troduced in this paper opens up the possibility of establishing similar results 
for other models where: (i) process-level diffusion approximations have been 
established; (ii) the stationary distribution of both the underlying queueing 
model and the associated diffusion limit is known to exist; and no rigorous 
link between (i) and (ii) is known. Notable examples of such systems are 
feedforward (acyclic) queueing networks and certain instances of multiclass 
networks for which stability conditions are known and Brownian approxima- 
tions have been rigorously established. One difficulty in directly applying the 
methods of this paper is the apparent lack of strong solution and Lipschitz 
continuity in Skorohod mappings corresponding to these models. Proving an 
interchange of limits for these classes of queueing networks is an interesting 
direction of future research that may benefit from the results we develop in 
the current paper. 
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APPENDIX: PROOFS 
A.l. Proofs of the main results. 
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Proof of Theorem 5. By definition of a geometric Lyapunov func- 
tion (23), we have that, for every x £ X, 

E x m~{t ))} < max{ 7 $(x),<Kto)$(z)l{$(z) < K}} 

<^(x)+(f>(t )K. 

Fix fceN and put <S>k( x ) := A k. Then, if < k 

$ k (x) - EspfcCEfo))] = - E x [$ k (E(t ))\ 

><S>(x)-E x M~(t ))} 

> (l-i)$(x)-</>(to)K 

> -Hto)K. 
On the other hand, if &k(x) > k, then 

$ fc (x) - E x [$ fc (H(t ))] > fc - E x [^ k (S(t ))] > 0. 

Thus, — E x [$jfc(S(io))] is bounded from below by —cf>(to)K for all 

x £ X . It is also bounded from above by k. Now, the monotone conver- 
gence theorem implies that E x $k(E(to)) j Ea.$(S(io)) as A; — >• co. So, $fc(x) — 
E x [$ fc (H(t ))] ->• $(x) -E x [$(E(t ))] ask^oo, for all xG/t. Since - 
E x [$jfc(H(to))] bounded below, we can appeal to Fatou's lemma to conclude 
that 



(37) 



($(x)-E x [$(E(t ))]Mdx) 



< liminf / (* fc (z) - E a .[$ fc (E(* ))]Mda:) = 0, 

fc— >0O J 



where the equality follows from stationarity of it. On the other hand, by 
definition of the Lyapunov function, we have that 

- E x [<S>(E(t ))} > (1 - 7 )$(x) - ^(to)^. 

Integrating both sides with respect to ir and using (37), we conclude that 

E^m) < f^, 

1-7 

which concludes the proof. □ 

PROOF of Theorem 6. Let 9 > satisfy (28). Introduce the function 
T : X 2 x R + -> R+, defined as 

T(xi,x 2 ,y) :=exp(y($(x 2 ) - $(sci)))- 
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Using a second-order Taylor's expansion around y = 0, we have, for some 
< ff < 0, 

T(si,x 2 ,0) = l + 0(*(a; 2 )-*(a:i)) 

+ y^fo) - $(^i)) 2 exp(0 , ($(x 2 ) - $(si))), 
<l + 0($(x 2 )-$(xi)) 

+ y ($(s 2 ) - $(zi)) 2 exp(0($(x 2 ) - $(xi)) + ). 

We now fix arbitrary x such that 3>(x) > if, set xi = x,x 2 = E(io) and take 
the expectation of both sides above to obtain 

E x [T(x,E(t ),e)} 

<l + 6M x mE(t ))-$(x)} 
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+ Y E x mZ(t )) - *(x))' exp(0(*(H(t o )) - Hx)) + )] 

to 2 , 

<i- 1 e + -L 2 {e,t ) 

(b) 

< 1 - 70/2, 

where (a) follows from (22) and the definition of L 2 , and (b) follows from 
(28). Note that the condition &{x) > K is equivalent to exp(0$(x)) > exp(0i^). 
Thus, exp(0<I>(x)) is a geometric Lyapunov function with parameters 70/2, 
io, exp(0i^). This concludes the proof of the first part of Theorem 6. 

To prove the second part, observe that, for the constructed Lyapunov 
function, 

(38) 0(i o )=su Pa; E a: [exp(0($(H(to))-$(x)))]=Li(0,fo). 

We use Markov's inequality, (38) and the bound (27) of Theorem 5, which 
yields 

P w (*(2(0)) >s) = P w (exp(0$(~(O))) > exp(0s)) 

< exp(-0s)E^[exp(0$(~(O)))] 

^ L 1 (0,t o )exp(^) 

< exp(-vs) —, 

- PV ' 1-70/2 

= (1 - 7 0/2)- 1 L 1 (0, t ) exp(-0(s - #))• 
This completes the proof of the second part of the theorem. □ 

Proof of Proposition 2. The key to the proof is the following result. 
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Lemma A.l. Let {X n (t) :t > 0) be the net input process for the GJN, 
E n , and let (x™(t) :t>0) denote its fluid counterpart starting at state z [i.e., 
x™(0) = z\. Then, there exists 9\ > such that 



(39) limsup sup n */ 2 E 



(40) 



and 



(41) 



limsup sup E sup 

n^oo ( z ,a,v)eX L(Kt<n 



sup \\X n (t)-x n z (t)\\\Q n (0) = (z,a,v) 

0<t<n 



< oo, 



\Q n (0) = (z,a,v) 



< oo 



limsup sup n 1 E 



sup \\X n (t)-x n z (t)\\'< 

0<t<n 



xe^(n _1/2 0i||Jr , (t)-a£(t)||) 



\Q n (0) = (z,a,v) 



< oo. 



The proof of this lemma is given in Section A. 2. 

From (8), (9) and (10), it follows that (Y n (■) , Q n (■)) jointly solve the 
Skorohod problem for the unreflected net input process X n (-). That is, 
Y n (-) = V 1 (X n (-)),Q n (-) = y 2 (X n (-)). Let (y n ,q n ) = {^ 1 (x n z )^ 2 (x n z )) be 
the solution of the Skorohod problem for the linear function x z (-). Fix 9 > 
and to > 0. Then, 

n^ 1 / 2 6(w'Q n (t n)-w'z) + 

< n- 1 / 2 9{w'Q n {t n) - w'q n (t n)) + + n~ 1 / 2 6(w'q n (t n) - w' z) + 

(42) <n -1 / 2 sup \w'Q n {t)-w'q n (t)\ 

0<t<t n 

(b) 

< C^r 1 / 2 sup \\Q n (t)-q n (t)\\ 

0<t<t n 

(c) 

^Rxddn- 1 / 2 sup \\X n {t) -x n z {t)\\. 

0<t<t n 

In the above (a) follows from from Proposition 1, from which we have that, 
for every z G Ri, 



(43) 



w'q n (ton) < w'z — min /x,(l — Pj)ton ) 
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w'z — min u~tr\K*\/n 
1<J<J J 

Thus, the second term on the RHS of the first inequality above is identically 
zero. The inequality (b) follows from setting C\ = maxwj, and (c) follows 
from (13) of Theorem 1. Thus, setting 9q = 6\/{R\Ci\JtQ ), where 9\ is iden- 
tified in Lemma A.l, concludes the proof of (31) and (32) for arbitrary choice 
of t . 

We now concentrate on proving (30). From (42) and (39) of Lemma A.l, 
we have that there exists a constant C% > such that 

E[\w'Q n (t n)-w'q n (t n)\\Q n (0) = (z,a,v)] 

<e[ sup \w'Q n (t)-w'q n {t)\\Q n (0) = {z,a,v) 

,0<t<t n 

holds for all (z,a,v) £ X, to > and all sufficiently large n, where we have 
substituted ton for n in (39). Combining the above with (43), we obtain 
that, for every to > 0, (z,a,v) £ X, and all sufficiently large n, 

E[w'Q n (t n)\Q n (0) = (z,a,v)} 

< max< C2\/ton, C2\/^o n + w'z — min HjtoK* \/n> 

= w'z + max< C2\/Uyn — w'z, C-iyfUsn — min UjtoK* 

We now set to as a function of C2,k* and min.j{//j}, so that C2\/to — 
mini<j< j fMjtoK* < —1. This is always feasible for large enough value of to 
since the expression above is quadratic in \JIq with negative coefficient in 
front of (y/to ) 2 . For said to, we obtain C2\/ton — mini<j< j fj,jtoK*y/n < —y/n. 
Then we set Co = C2y/to + 1- Consequently, for values of z such that w' z > 
CQ\/n, we have that C2\Ao n — w'z < —\fn. We conclude that for our con- 
stants to , Co and C2 , we have that 

E[u/Q n (t n)|Q n (0) = (z,a,v)] <w'z-y/n~ 

whenever w'z > CQ\fn. This proves (30). □ 

Proof of Theorem 7. We apply Proposition 3. Fix a constant < 
c < 1. The proposition implies that, for 6 n = cOqu -1 ^ 2 , we have 

(44) limsupLi(0 n , nto) < 00 
and 

(45) limsupn _1 L2(^n, nto) < cxd. 
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This implies that 

lim.supn~ 1 / 2 9 n L2(0 n ,nto) = limsupc9on~ 1 L2(9 n ,ntcj) < 1, 

n— >oo n^oo 

provided that c is sufficiently small. This means that condition (28) of The- 
orem 6 holds for all sufficiently large n. 

Applying Theorem 6 and using n 1 / 2 9 n = c9q, we obtain, for every s > 0, 

P,„KQ n (0)/^ > s) < L ^\ nto J eM-0n(sn^ - con 1 / 2 )) 

1 — cOq/2 

= l-c9 /2 eM-c9 (s-c )). 

From (44), there exists a sufficiently large no and a constant C3 > such 
that Li(9 n , nto) < C3 for all n > uq. Then for all s > cq and all n > uq, 

V* n (w'Q n (0)/^ >s)< °\ exp(-c9 (s - c )). 

1 - c9 /2 

Setting C\ and c\ appropriately, we obtain the result. □ 

Proof of Theorem 8. Corollary 1 implies that, for every subsequence 
7r nk , there exists a weakly convergent subsequence Tr nk i . Let ir* denote any 
weak limit of this subsequence. We claim that every such limit measure ir* 
is equal to 7Trbm> implying that all such weak limits are the same. It then 
follows that 7r n => ttrbm as n — > oo. Thus, we are left to show that every 
converging sequence 7r nfe * ir* converges to 7Trbm . For notational simplicity, 
we denote n nk i simply by n n . 

We apply Theorem 4. Fix an arbitrary t > 0. Since 7r n converges weakly 
to tv*, then the process (Q n {nt')/y/n: < t' < t) conditioned on Q n (0) ~ vr n 
converges weakly to a (/?, T,I-P')-KBM (Z(t') :0<t' <t) with initial distri- 
bution given by ir* . In particular, Q n (nt)/y/n converges weakly to Z(t) with 
Z(0) distributed according to tt* . But from stationarity of 7r n , Q n {nt)/\fri 
is distributed as 7r n , which, by assumption, converges to n* . Therefore, Z(t) 
is distributed as tt* and tt* must be the unique stationary distribution: 

7T*=7r R BM- D 

Proof of Corollary 2. It suffices to prove the asserted uniform inte- 
grability. Fix 5 > 0. By Theorem 7, we have that, for all sufficiently large n, 



E[exp(n- 1 / 2 5w'Q n (0))}= / F(n" 1/2 w'Q n (0) > <T X logx) dx 

Jo 

< / exp(— ciS^ 1 logx + C2) dx < 00, 
Jo 

provided that c\5~ l > 1. Setting 9 = 8w, we obtain the result. □ 
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Proof of Theorem 9. The proof repeats the argument used in [33]. 
It is first shown in [33] that, when Q n (0) = 0, 



(46) lim E 



sup Df >h (t')/n 

0<t'<nt 



0. 



for every j G J . This convergence is shown by bounding the supremum by 
a certain linear functional of sup <t/< n £ ||Q n (t')||/n. Then a random time 
change theorem is applied to show that the distribution of sojourn times 
in station j at time nt' becomes "indistinguishable" from the distribution 
of the sum of the workloads in every individual station j' multiplied by hi . 
The latter converges weakly to hjifij>Zj(t') by Theorem 4, and the result fol- 
lows. But (46) also holds in our setting when Q n (0) is not zero, but instead 
Q n (0)/y/n =>• vr^. Indeed, by Theorem 4, (Q n (t') / y/n: < t' < nt) converges 
weakly to a corresponding RBM with Z(0) ~ ttq. Then sup 0<t , <nt \\Q n (t')\\/y/n 
converges weakly to sup 0<t / <t ||Z(i')||, with initial distribution 7Tq. This im- 
plies that, for every e > 0, 

lim pf sup \\Q n (t')\\/n>e] = lim pf sup ||Q n (t')||/v^ > ey/n) = 0, 

n^oo Vo<t'<ni / n ^°° \0<t'<nt / 

since E^* [sup < i /<t ||Z(i')||] < oo. Repeating the argument in [33], (46) then 
holds as well and the statement in the theorem follows. □ 

A.2. Proofs of auxiliary results. 

Proof of Lemma A.l. We first show that (40) implies (39) and (41). 
Given 9\ > 0, select C large enough so that exp(^ix) > x 2 > x for all x > C. 
Replacing ra" 1 / 2 sup 0<t<n \X n (t) — x™(t)\ with max(C, n" 1 / 2 sup 0<t<n ||A n (t) - 
rr™(t)||), we obtain that (39) and (41) hold, provided that 



lim sup sup E 



sup exp(2^max(C7,n- 1 / 2 ||X n (t) - x?{t)\\)) 

0<t<n 



\Q n (0) = (z,a,v) 



< oo. 



But this bound holds if and only if (40) holds. Thus, it suffices to prove (40). 

The proof of (40) relies on considering tail bounds on the uniform devi- 
ations of X n from Xg. We expand as (aj"i(t), . . . ,x™ j{t)). Note that, 
because X n (0) = a£(0) = Q n (0) = z, we have, for all j = 1,'. . . , J, 

X?(t)-x?j{t) = (A?{t)-ap) 

+ y: Pimm))-^m))-{sAm))-M B W)) 

l<i<J 
Ki<J 
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where we recall that the process S(t) is independent from n. 

Fix j, 1 < j < J. We then have, for all integer n > 1 and (z, a, v) 6 A", 

P( sup exp(^ 1 / 2 |X ? n (t)-x"-(t)|)>u|Q n (0) = (z,a,?;) 

\0<t<n 

= p( sup \X?(t)-x%At)\ >6- 1 n 1/2 logu\Q n (0) = {z,a,v) 

V 0<t<n 

<p( sup L4?(i)-a"i| >(8~ 1 n 1 / 2 logu)/3\Q n ( y 0) = (z,a,v)) 

\0<t<n / 

(47) +pf sup \Sj(m(t)) - (ijB?(t)\ >(^ 1 n 1 / 2 logu)/(3(J+l)) 

\0<t<n 

|Q n (0) = ( 2 f J a,i;; 

+ p( sup |^.(^(^(t)))-^S i (i?™(t))|>(^ 1 n 1 / 2 log U )/(3J) 

\0<t<ra 

|Q"(0) = (-W) 

for all it > 1. Since £?"(t) — £?™(s) < (t — s) for all t > s, it suffices to bound 
the uniform deviation of Sj(t) from [ijt, and Rj(Sj(t)) from PijSi(t) in the 
latter two terms on the right-hand side above, for all i,j = l,...,J. We 
now proceed to bound these terms using the so-called coupling probability 
bounds that are the main building blocks in strong approximations (cf. [8], 
Chapter 2.6). Introduce the following J(J + 2) constants: 

r °j= a j c l,j for j = l,...,J, 

(48) I'/'' I' A.j for j = l,...,J, 

Tfcj = Pfcj(l - Pkj) for k,j = l,...,J. 

We now use the following lemma, which provides bounds on the probability 
of deviation between the stochastic primitives and certain Brownian motion 
processes that are constructed on the same underlying probability space. 

Lemma A. 2. There exist J + 2 mutually independent standard J -dimen- 
sional Brownian motions W k = (W k (t) : t > 0) , k = 0, . . . , J+ 1, such that 

sup \Aj(t) - a 3 t - (lQ) 1/2 W}(t)\ > Cgj logn + u 

0<t<n 



|Q(0) = (z,0,0))<C 2 ,e- Co ^, 
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P( SUp \Sj(t) - fijt 
\0<t<n 

- (r/ +1 ) 1 /2 W /+i( t )| >C) +X j\ogn + u 
f( sup \B^(lt\)-p kj t-(Xkj) 1/2 W}(t)\ >C7^-Iogn + u 

\ 0<t<n 

|Q(0) = (*,0,0)) <Cl je - c ^ u , 

for j = 1, . . . , J and fe = 1, . . . , J, where C\ ■ , Cf ^ and Cfc j , k = 0, . . . , J + 1, 
j = 1, . . . , J, are /iraie positive constants independent of u and n. 

Note that the first identity is given for the limiting process A(t) and not 
the indexed process A n (t). But since A n (t) = A(a~ 1 a'jt),a'jt = (aj afflatjt 
and a™ < ctj, then sup < i < n \ A^(t) — Ojt\ < sup <t< n \ Aj(t) — atjt\ . Therefore, 
it suffices to obtain an appropriate bound on sup Q<t<n \Aj(t) — ajt\. Note 
also that in the first expression above conditioning on Q(0) = (z, 0,0) is 
equivalent to conditioning on ao = 0, since the arrival process is indepen- 
dent of the service time process and the queue length vector at time t = 0. 
Similarly, for the second expression, the conditioning is equivalent to condi- 
tioning on vq = 0; for the third expression, conditioning on Q(0) = (z,0,0) 
is redundant. 

Using Lemma A. 2, we can bound the three terms on the right-hand 
side of (47). Since the proof is essentially identical (with obvious modifi- 
cations) in all three cases, we will provide a complete argument only for 
the first term on the right-hand side of (47), that is, P(sup 0<t<n | Aj(t) — 
ctjt\ > (9~ l n l l 2 logn)/3|Q(0) = (z,a,v)), omitting the other cases. Let r = 
(# _1 n 1 / 2 logit)/3. Then, we have 

P( sup \Aj(t) - ajt\>T\Q(0) = (z,a,v) 

\0<t<n 

= P( sup \Aj(t) - ctjt\ > t|oj(0) = aj 

\0<t<n 

<P( sup \A j (t)-a j t\>T\a j (0)<T/(2a j ),a j (0) = a j 

\0<t<n 

(49) x P(a i (0) < r/(2o !j )|a j (0) = a,) 

+ P(a i (0) > 7-/(20,-) I (0) = a,j) 
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< 



'( sup \Aj(t) - ajt\ > r\aj(0) < r/(2a j ),a j (0) = aj ) 

\0<t<n J 



+ P(o i (0) > T/{2a j )\a j (0) = a j ). 

We now bound each of the two terms appearing in the sum above, starting 
with the second term. Using Markov's inequality, we have that 

P(o,(0)> T/(2a j )\a j (0) = a j ) < E[exp(0Oj(O))|a,,(O) = a,] exp(-0T/(2ay)) 

= E[exp(0a i (O))|a j (O) = aj]u~^~ lnl '\ 

using the definition of r. Recall that the distribution of Oj(0) conditioned on 
Oj(0) = aj is the same as the distribution of aj(l) — aj conditioned on aj(l) > 
aj, and from (1), we have that, for 9 < 6*, sup a gK+ E[exp(#Oj(0))|aj (0) = 

aj] < oo. When n > 2 2 (6aj) 2 , we have u -( 6a j)~ lnl/2 < u ~ 2 (the choice of 2 
as a lower end for the interval of integration below will become clear in what 
follows). This implies that, when n > 2 2 (6aj) 2 and 6 <9*, we have 



limsup sup / P(aj(0) > T/(2otj)\aj(0) = aj) du 

(50) < limsup sup E[exp{6a j (0))\a j (0)=a j ] [°° u ~^ a ^ lnl/2 du < oo. 

Now we estimate the first term in (49). For every r' < r/(2aj), condi- 
tioning on aj(0) = t' , we have Aj(0,t) = 1 + Aj(r' ,t). Observe that the dis- 
tribution of Aj(r' ,t) is the same as Aj(0,t — r') conditioned on having the 
residual interarrival time a, (0) = 0. This follows since, by assumption, there 
was an arrival at time (time r' before the shift). Thus, we obtain 

sup \Aj(t) — ctjt\ > r\aj(0) = t ,cij(0) = aj 

0<t<n 

sup \1 + Aj(t - t') - ctjt\ > r\aj(0) = 

r'<t<n 

<P( sup \A j (t-T / )-a j (t-T / )\ + l + a j T / >T\d j (0) = 

\T><t<n 

<P( sup \A j (t)-a j t\ + l + a j T'>T\a j (0) = 

\0<t<n 

<P( sup \A j (t)-a j t\>(T/2)-l\a j (0) = o), 

\0<t<n J 

where in the last inequality we have used the fact that ayr' < r/2. Since the 
bound above is uniform over r' in the range [0, r/{2aj)], it follows that 

P( sup \A j (t)-a j t\>T\a j {0)<T/2,d j {0) = a j ) 

\0<t<n J 
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<P sup \A j (t)-a j t\>{T/2)-l\a j (0)=0). 

\0<t<n J 

Now, recall that r = (9~ 1 n 1 ^ 2 logu)/3, thus, plugging into the above bound, 
we get 



0<t<n 

< _ 

1 0<t<n 



sup \Aj(t)-ajt\ > (0~V/ 2 logu)/6- 1|%(0) = 

sup \Aj(t) - ajt - (TQ)V 2 W?(t)\ > (e- 1 n 1 / 2 log«)/12 - 1 

3<t<n 

(51) + p(V 1/2 sup \{T^) 1/2 W^(t)\ ^(O-Hoguyu 



0<t<n 

< n c i exp(-(cj /9)y/n\ogu + 1) + Cj exp(-c^~ 2 (logn) 2 ), 



\2\ 



(52) = en^u-te/*)^ + Cj exp(-c^- 2 (log^ 

for some finite positive constants Cj,C'j and Cj,c'j, where the last step follows 
from the first displayed equation in Lemma A. 2 and standard tail bounds 
on the maximum of a Brownian motion (see, e.g., [8], Lemma 1.6.1). For the 
second term in the RHS of (52), we have 



oc 



2 



Let no be the smallest integer such that (cj/9)y/no> 1. Then, for the first 
term on the RHS of (52), we have 

/•OO /"OO 

limsup / n c iu~ {c 3 /d)V7l du< sup / n c %^/ e )^"cfri 

n— >oo J 2 n>nr,J2 



= SUP ; ,-, _ : < OO. 

n>no (( Cj /^)^- l)2fe/ e) ^"~ 1 

By the finiteness of the two bounds above, we obtain 

limsup/ P( sup -oyi| > (0 -1 n 1/2 logit)/6 - 1 

n-*oo J 2 \0<t<n 

(53) 1^(0) < r/(2a i ),a j (0) = a,- J du 

< oo. 

Putting this together with (49) and (50), we have 



lim sup sup E 

1-+00 ( 2)0 ,v)eA' L0<t<n 



sup exp(tm 



- 1 / 2 |A i (t)-a i t|)|Q(0) = (z J a,i;) 
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roc / 

<2 + limsup sup / PI sup \Aj(t) - ctjt\ > {6- l n l l 2 \ogu)/?> 

n (z,a,v)eXJ2 \0<t<n 

(54) 



\Q(0) = (z,a,v)j du 

< oo. 

Repeating these arguments for the service completion processes, we obtain 
similar bounds for the second term on the RHS of (47). To obtain the bound 
on the third term, 

pf sup \R)(Si(t)) - Pij Si(t)\ >(6- 1 n 1 / 2 logu)/(3J)\Q n (0) = (z,a,v) 

\0<t<n 

we set again r = # _1 n 1//2 logu/3 and note that this term is bounded by 
Pf sup \R)([t\) - Pij t\ >T/J\Q n (0) = (z,a lV )) 

\0<t<n+r/(J+l) / 

+ pf sup \Si(t)- mt\>T/(J + l)\Q n (0) = (z,a,v)). 

\0<t<n J 



(55) 



We have shown above that, for the second probability in the sum above, the 
corresponding event sup 0<i<n \Si(t) — /ijt| > r/(J + 1) satisfies the bound 
of the form (54) with Si and fa replacing Aj and ctj, respectively. We now 
bound the first term in the sum above. Note that the corresponding event is 
independent from the event Q n (0) = (z, a, v) since it only involves the routing 
process. We apply Lemma A. 2 and, in particular, consider the corresponding 
coupled Brownian motion W(t) with variance term r,j =pij(l —pij). We 
have 



pf sup \R)([t\)- Pl] t\>r/j) 

\0<t<n+r/(J+l) / 

<pf sup |m(Ltj)_p.^-(r M -) 1/2 Ty(n + r/(J + l))|>r/J 

\0<t<n+T/(J+l) 

+ Pf sup |(r M -) 1/2 ^(n + r/( J + 1))| >r/j). 

\0<t<n+T/(J+l) J 

To bound the first term, we use the last part of Lemma A. 2 and rewrite the 
term as 

sup \R)([t\ ) - Pij t - (TijWwin + r/( J + 1))| 

0<t<n+r/{J+l) 

> Cl d log(n + r/( J + 1)) + (t/J - C\ 4 log(n + r/( J + 1))) 
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< Cfj exp(-c itj T/J + CijCjj log(n + r/( J + 1))) 

= (n + c^n 1 / 2 logn) C2 exp(— c^ra 1 / 2 logu), 

where the last term is obtained by expanding r as _1 n 1 / 2 logu/3 and in- 
troducing the necessary constants c' k ,k = 1,2,3. For all sufficiently large n, 
the last term is integrable with respect to u and with the boundaries 
uniformly in n. Namely, 

limsup / (n + c^n 1 / 2 log u)° 2 exp(— c^n 1 / 2 log u) du < oo. 

n—>oo J 2 

To bound the second probability in (55), we again use standard tail bounds 
on the maximum of a Brownian motion as in (51). The derivation of the 
bound is very similar, the only difference being that the upper limit n in the 
supremum expression is replaced byn + r/(J+l) = n + 6~ l n l l 2 logu/ (3( J + 
1)). An easy check shows that the expression is still uniformly integrable in 
u and the bound similar to (54) is obtained. We conclude that 



limsup sup / P( sup \R l j{Si{t)) — pijSi{t)\ 

n-'oo ( Z nv)£XJ2 \0<t<n 



> (0-V /2 logu)/(3J)|Q n (O) = (z,a,vfj du < oo. 
To summarize, for any < 6 < 9* , we have 

9n~ 1 / 2 \X^(t) - x](t)\)\Q n {0) = (z,a,v) < oo. 

Setting 6*i < 8*/ J, we obtain 

^ 1 ' 2 \\X n {t) - x n {t)\\)\Q n {<d) = (z,a,v) < oo, 
which concludes the proof of (40). □ 



lim sup sup E 

n-*oo ( z ,a,vjeX 



lim sup sup E 



Proof of Lemma A. 2. The results of this lemma are immediate appli- 
cations of functional strong approximations for random walks and the asso- 
ciated renewal processes, [7, 24, 25] (see also [11], Chapter 5). In particular, 
the following approximation results are established in the above referenced 
papers. 

Theorem 13. Let £j, i = 1,2, . . . , be an i.i.d. sequence such that E[e <5 ^ 1 ] < oo 
for all sufficiently small 5 > 0. For every t G define X(t) = J2i<[t\ £i 
and the associated renewal process Y(t) = max{/c : J2i<i<k < There ex- 
ist standard Brownian motions W\(t), Wzft) defined on the same probability 
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space as £j, and constants C±, C2, C3 > suc/i i/iai ; /or every T > 1, it > 0, 

(56) p( sup |X(t)- l /- 1 t-o-VF(t)|>CilogT + n) <C 2 e- C3n , 

\0<i<T / 

(57) Pf sup \Y(t) -vt- {vfl 2 oW(t)\ >CilogT + u) < C 2 e~ C3U , 

\0<t<T J 

where v = 1/E[£i] and a is the standard deviation o/£i. 

The results of the lemma follow immediately from this theorem. The 
bounds corresponding to the arrival processes Aj(t) and service processes 
Sj(t) follow from (57) directly by setting v = ay and v = respectively, 
and observing that v 3 / 2 a = otjC 2 .^ for the jih arrival process and = fijc 2 • for 
the jth service process. The last bound corresponding to the routing process 
Rj([t\) follows from (56) by setting v = Pkj and a = (T^j) 1 ' 2 . □ 
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